Models of electrolyte solutions from molecular descriptions: The example of NaCl 

solutions 
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We present a method to derive imphcit solvent models of electrolyte solutions from all-atom 
descriptions; providing analytical expressions of the thermodynamic and structural properties of 
the ions consistent with the underlying explicit solvent representation. Effective potentials between 
ions in solution are calculated to perform perturbation theory calculations, in order to derive the 
best possible description in terms of charged hard spheres. Applying this method to NaCl solutions 
yields excellent agreement with the all-atom model, provided ion association is taken into account. 
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Since the pioneering works of Debye, Hiickel, and 
Onsager, electrolyte solutions have been commonly 
described by continuous solvent models, for which 
the McMillan-Mayer theory [ij provides a rigorous 
statistical-mechanical foundation. Within that level of 
description, simple phenomenological models such as the 
primitive model (PM), for which the ions are assimi- 
lated to charged hard spheres 0, can lead to explicit 
formulas for the thermodynamic and structural proper- 
ties (e.g., with the help of the mean spherical approxima- 
tion (MSA) or the binding MSA (BIMSA) [|). These 
models are the most practical to use Q , since they allow 
for a direct link between the experimental measurements 
and the microscopic parameters of the system. Never- 
theless, they ignore the molecular structure of the sol- 
vent. Consequently, they cannot properly account for 
the complex specific effects of the ions, which appear in 
numerous biological, chemical, and physical interfacial 
phenomena 0,(3 j without further developments. 

An alternative procedure consists in carrying out 
molecular simulations, where both the solvent and solute 
are treated explicitly. After a rigorous averaging over 
the solvent configurations, a coarse-grained description 
of the ions, which still includes the effect of the solvent 
structure, can be obtained f8l-[T]|. However, this set of 
methods is purely numeric; they do not provide any an- 
alytical expression for thermodynamic quantities. They 
arc therefore restricted to simple geometries [H,[ll] (bulk 
solutions or planar interfaces). The description of com- 
plex systems, such as porous or electrochemical materi- 
als, is still based on continuous solvent models 

In this letter we present a method aimed at bridging 
the gap between analytical and numerical approaches. It 
is based on the application of liquid perturbation theory 
(LPT) [l^ to effective ion-ion potentials extracted from 
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molecular dynamics (MD) results. Different approxima- 
tions of the PM are employed for the case of NaCl elec- 
trolyte solutions: a two component model (MSA2), that 
only takes free ions into account, and two different three 
component models (MS A3 and BIMSA3), which include 
a third species (the contact ion pair). As we proceed 
to show, LPT allows us to select the best simple model 
which accurately accounts for the thermodynamics and 
the physical-chemistry of the system. 

The first stage consists in calculating the McMillan- 
Mayer effective ion- ion interaction potentials X^j^(r), by 
inverting the radial distribution functions (RDF) gij{r) 
obtained by MD. The simulations were carried out on 
a box of 2000 water molecules and 48 NaCl pairs us- 
ing the same interaction potentials as in reference [19. 
This setup corresponds to a concentration of 0.64 mol P . 
NPT ensemble sampling at standard pressure and tem- 
perature was enforced, with a time step of 1 fs and a 
pressure bath coupling constant of 1 ps. An equilibration 
run of 0.25 ns was followed by a production run of 0.6 ns 
for five different initial configurations. The averages of 
the resulting RDF were then used for the potential inver- 
sion via the HNC closure [l5[ . These effective potentials 
are assumed to be concentration independent and will be 
used for simulations at all concentrations. 

Subtracting the long-range Coulombic potential 
Vlj^'{r) (which depends on the dielectric constant of the 
solvent) from Vfj^{r), we obtain the short-range contri- 
bution V^^{r) to the effective potentials. These are given 
in Fig. [1] (species 1 and 2 refer to Na"*" and Cl~ free ions, 
respectively). All the short-range potentials exhibit os- 
cillations corresponding to the solvent layering between 
the ions, but this effect is particularly important for the 
cation-anion interaction: a considerable potential barrier 
(> 2k'QT) separates the first two attractive wells. To 
serve as a reference, Monte Carlo (MC) simulations were 
performed with these effective potentials; a comparison 
between MD and MC RDF is also provided in Fig.[TJ The 
excellent agreement between both sets of RDF validates 
the HNC inversion procedure |17| , and allows us to com- 
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FIG. 1: Effective McMillan-Mayer short-range pair potentials 
extracted from explicit solvent simulations using the HNC 
closure, (a) Cation anion, (b) cation cation, (c) anion anion, 
(d) cation anion RDF obtained from explicit solvent MD and 
implicit solvent MC simulations. 

pute all ion thermodynamic properties through implicit 
solvent MC simulations. 

The second stage of our coarse-graining procedure con- 
sists in applying LPT, in order to deduce the best ana- 
lytical model of electrolyte solutions which reproduces 
this molecular description. The principle of LPT is to 
describe the properties of a given system in terms of 
those of a well known reference system, with the differ- 
ence between them treated as a perturbation in the ref- 
erence potential. Assuming pairwise additive potentials, 
Vij — Vjj^^ + AVij , a first-order truncated expression for 
the free energy density of the system is obtained, 

Pfv < /3/i°) + E P^PJ J dr 9^f{r)AV,, (r) (1) 

which depends only on the free-energy density fv^^ and 
RDF of the reference fluid, with ^ = (fceT)-! and 
Pi the concentration of species i. The Gibbs-Bogoliubov 
inequality [HI ensures that the right-hand side of Eq. ^ 
is actually a strict upper bound. Once a reference system 
has been chosen, the expression on the right-hand side of 
Eq. ([1]) must be minimized with respect to the parameters 
defining the reference. This procedure yields the best 
first-order approximation to the free energy of the system 
under consideration. 

For a system of charged particles in solution, the nat- 
ural reference is the PM, defined in terms of the charge 
and diameter (ai) of each species. In this case, the per- 
turbing potentials are just the short-range effective po- 
tentials computed above (Al/^j = V^^^). We use the 
MSA 0] solution to the PM, since it provides analyti- 
cal expressions for both the free energy and the RDF. 
The perturbation term is evaluated using an exponential 
approximation to the RDF obtained within the MSA, 
g{r) = exp [gMSA{r) — 1], which removes any unphysical 
negative regions and improves the comparison with HNC 
calculations. 



FIG. 2: (Color online) (a) Osmotic coefficient "1? in the 
McMillan- Mayer frame of reference, (diamond) MC simula- 
tions, (dot dashed) MSA2, (dot) Debye Hiickel Limiting law 
(DHLL), (cross) experiments (Ref. [ig] with the McMillan- 
Mayer to Lewis Randall conversion), (b) Minimization diam- 
eters, (dot dashed) MSA2 and (diamond) MSA-fit. 

We first used LPT for a two-component system (Na+ 
and CI" free ions) within the MSA (model MSA2), for 
concentrations ranging from 0.1 to 2.0 moll~^. The mini- 
mization leads to almost constant diameters on the whole 
range of concentration: cti = 3.67 A and (T2 = 4.78 A. 
As shown in Fig. [2l these parameters yield osmotic co- 
efficients close to MC calculations only at very low con- 
centration, i.e., c < G.lmolP^ (experimental values are 
given for indicative purposes only, since a perfect model 
will exactly match the MC results). For molar solutions, 
the LPT results differ considerably from MC calculations. 
This discrepancy can easily be understood by comparing 
the diameters found within the MSA2 calculation with 
the effective potentials given in Fig.[TJ The anion/cation 
contact distance obtained within the MSA2 calculation 
is 4.2 A, which is in the region of the second minimum of 
the effective potential and corresponds to the situation 
where there is a single layer of water molecules between 
the ions. The first minimum of the potential, which cor- 
responds to the contact ion pair (CIP) is thus completely 
ignored by the MSA2 calculation. If the MSA diameters 
are directly fitted to reproduce the MC osmotic pres- 
sure, much smaller values are obtained. These MSA-fit 
hydrated diameters, which are compared to the MSA2 
diameters in the bottom part of Fig. [2l are averages of 
the CIP and the solvent-separated ion pair. 

To overcome this difficulty, we have explicitly intro- 
duced the CIP in our model (species 3). Straightforward 
calculations, based on a characteristic-function formal- 
ism, allow us to define an equivalent model in which 
the free ions and the CIP are explicitly taken into ac- 
count [l^, [lOl- We apply this formalism by defining a 
pair as an anion and a cation at a distance less than 
4 A, which corresponds to the position of the effective 
potential maximum. The interaction between free, like 
charges in this new system remains unchanged, and the 
cation-anion interactions are easily approximated by ex- 
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FIG. 3: EfTective pair potentials derived for MSA3 and 
BIMSA3. (a) Cation anion (dashed line: without taking the 
pair into account), (b) pair cation, (c) pair anion, and (d) pair 
pair. The internal potential of the pair l3Vint{r) is set equal 
to I^Vif^lr) for distances less than 4 A. 



trapolating the original potential at the barrier separat- 
ing pairs from free ions (as shown in Fig. [3|). We assume 
that the interaction potential is averaged over the rota- 
tional degrees of freedom of the CIP and thus pairwise 
additive. Hereafter, the quantities referring to such a 
three-component model are written with a tilda symbol. 
The short-range potentials involving the pair can be de- 
rived, in the infinite dilution limit, from an average of 
the contributing ion interactions. In Fourier space, 

Vi^ik)^wik/2)[V,^^+V,T]ik), ^ = l,2 (2a) 
Vi^'^ik) = w{kl2f [ySR + ySR ^ 21/SR] (fc) (2b) 

where w{r) is the pair probability distribution 

w{r) = A'-^e-^^-'^'') (2c) 

Vint(?') is the internal part of the pair potential (see 
Fig. [3|), and Kq is the association constant, defined as: 

/>oo _ 

A'o = / dr47rr2e-'^^-'(''' = 0.43 L.mor^ (3) 
Jo 

The excess free-energy density of the original system 
is that of the three component mixture plus a 
correction term 

/3/r = /3/r-P3lni^o, (4) 

which is due to the change in standard chemical potential 
between the two component and three component mod- 
els. It should be noted that the fraction of pairs is now an 
additional parameter in the minimization scheme, which 
serves to ensure chemical equilibrium. Within this rep- 
resentation, the pair can be modeled as a hard sphere 
(MS A3) or as a dumbbcU-like CIP (BIMSA3) 0. Since 



FIG. 4: (Color online) Excess free-energy density Pf^^ as 
a function of the square root of the concentration ^/c. (dia- 
mond) MC simulations, (dot dashed) MSA2, (dashed) MSA3, 
(solid) BIMSA3, (dot) DHLL, and (cross) experiments. The 
inset gives the fraction of pairs (MSA3, BIMSA3) as a func- 
tion of ^/c. 

we have no additional information, we consider only sym- 
metric dumbbells. Furthermore, since analytic expres- 
sions for the RDF within BIMSA are not known, we ap- 
proximate the dumbbell as a hard sphere when comput- 
ing the perturbation term (this is not necessary for the 
reference term, since an expression for the free energy 
is available). Let CTc be the diameter of the cation (an- 
ion) within the dumbbell, the diameter of the hard sphere 
representing this dumbbell is taken to be (T3 = ^-^g^[2l|. 

Using these two reference systems, the three- 
component MS A3 and BIMSA3, we obtain results in 
much better agreement with the MC simulations, as 
shown in Fig. |4l The diameters obtained for species 1, 
2, and 3 are 3.65, 4.79, and 5.76 A for MSA3 and 3.69, 
4.75 and 6.19 A for BIMSA3. The free ion diameters are 
similar for MSA2, MSA3, and BIMSA3. The pair diam- 
eter is smaller when modeled as a hard sphere (MSA3) 
than when modeled dumbbell (BIMSA3). At high 
concentration (about 1 molP^), the MSA3 overestimates 
the free energy, because the excluded volume repulsion 
becomes too important for the pairs to be represented as 
hard spheres. The BIMSA3 model is the closest to the 
MC simulation results. It is worth noting that even at 
the lowest concentration considered, the fraction of pairs 
(shown in the insert of Fig. |4]), although less then 5%, 
has a non-negligible effect on the thermodynamics of the 
system. 

This procedure also provides an accurate description of 
the structure over the whole range of concentrations. A 
development similar to the one that leads to Eq. ^ de- 
rives the average unpaired RDF from the corresponding 
paired quantities: 

PiP]9i]{k) = P3w{k) (1 - 6ij) + p.iPjgij{k) 

+ p3w{k/2) [pigsi + pjgsj] (fe) (5) 
+ p^[w{k/2)fg33{k) 
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FIG. 5: (Color online) RDF obtained from MC simulations 
(diamond), BIMSA3 (solid line), and MS A- fit (dot dashed) 
at two concentrations. 

The RDF obtained within BIMSA3 arc compared with 
the MC and MSA-fit results in Fig. [3 Our BIMSA3 
model accounts for the strong molecular peak of the CIP 
and provides the correct distances of minimal approach; 
whereas the naive MSA-fit procedure ignores the former 
and gives poor estimates for the latter. At larger sep- 
arations, the BIMSA3 results do not reproduce the os- 
cillations observed in the MC simulations, but the cor- 
responding energy oscillations in the effective potentials 
are less than fceT. In addition, the perturbation term 



of the BIMSA3 appears to be negligible compared to the 
reference term for concentrations less than 1 moll~^. The 
perturbation can then be omitted to obtain a fully ana- 
lytical theory, determined by the hard sphere diameters 
and the pair fraction given by LPT; with the free energy 
and the RDF given in terms of the BIMSA and MSA so- 
lutions, as described above. While the procedure we have 
followed uses two different approximations for the refer- 
ence and perturbation terms (MSA vs BIMSA), these are 
known to be accurate for the systems under consideration 
and do not appear to be inconsistent with each other. 

To conclude, we have combined MD simulations with 
LPT to construct simple models of electrolyte solutions 
which account for the molecular nature of the solvent. 
The final result is fully analytical and it yields the ther- 
modynamic and structural properties of the solution, in 
agreement with the original molecular description. The 
methodology can in principle be adapted to any molecu- 
lar description of the system (MD simulations involving 
interaction potentials accounting for polarization effects 
or Car-Parrinello MD simulations for example) as long 
as the ion-ion RDF are known. It can also be generalized 
to study interfaces. The method appears to be a promis- 
ing approach toward the description of the specific effects 
of ions, especially for complex systems whose modeling 
requires an analytic solution. 

The authors are particularly grateful to Werner Kunz 
for fruitful discussions. 
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